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In the exact Kohn-Sham density- functional theory (DFT), the total energy versus the number of 
electrons is a series of linear segments between integer points. However, commonly used approximate 
density functionals produce total energies that do not exhibit this piecewise-linear behavior. As a 
result, the ionization potential theorem, equating the highest occupied eigenvalue with the ionization 
potential, is grossly disobeyed. Here, we show that, contrary to conventional wisdom, most of 
the required piecewise-linearity of an arbitrary approximate density functional can be restored by 
careful consideration of the ensemble generalization of DFT. Furthermore, the resulting formulation 
introduces the desired derivative discontinuity to the exchange-correlation potential. All this is 
achieved while neither introducing empiricism nor changing the underlying functional form. The 
power of the approach is demonstrated on benchmark systems using the local density approximation 
as an illustrative example. 
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Density functional theory (DFT) is an approach to the 
many-electron problem, which is based on mapping the 
interacting electron system into a non-interacting one. 
Exact in principle, the mapping depends in practice on 
an approximate expression for the exchange-correlation 
(xc) energy as a functional of the density, E xc [n{r)\ [1— 
5]. While the exact form of E xc [n(r)] remains unknown, 
many constraints it has to satisfy have been formulated. 
Of particular interest here is the piecewise-linearity prop- 
erty: Using a zero-temperature ensemble of integer elec- 
tron states [6, 7], the realm of DFT has been extended to 
fractional electron numbers (N = Nq + a, where N e N 
and a G [0,1]). It has been shown [8] that the total 
ground-state energy, E, is given by 

E(N) = (l-a)E(N )+aE(N + l). (1) 

An important practical manifestation of the formal 
piecewise-linearity condition [8-12] is the relation be- 
tween the highest occupied orbital energy, e^ OJ and the 
ionization potential, / = E(No) — E(No + 1). If piccewise 
linearity is maintained, it can be shown that Eho = —I, 
a result known as the ionization potential (IP) theo- 
rem [8, 13]. 

The above discussion establishes the importance of 
approximating the xc-functional such that piecewise- 
linearity in the total energy is obtained. However, it has 
long been known that commonly used functional classes, 
such as the local density approximation (LDA), the gen- 
eralized gradient approximation (GGA), or the conven- 
tional hybrid functional approximation, grossly disobey 
this condition in their standard implementation. Instead, 
a typically convex E(N) curve is obtained (see, e.g., 
[9, 12, 14-18]) and, correspondingly, the discrepancy be- 
tween Sho and —I can easily be as large as a factor of 
two [19]. 

Two main approaches have emerged in response to this 
problem. In one approach, various correction terms are 
imposed on existing underlying xc- functionals [20-27] . In 



another, piecewise linearity is explicitly enforced in the 
construction of novel range-separated hybrid function- 
als [28-33]. 

The above considerations on piecewise-linearity, or lack 
thereof, are all based on a description of fractional- 
electron systems by insertion of a density n(r), leading 
to a fractional total number of electrons, into a density 
functional developed originally for pure states. One may 
question whether this straightforward application is at all 
optimal. Indeed, Gidopoulos et al. [34] have observed, in 
the context of an excited-state ensemble, that straight- 
forward application of the Hartree term leads to an un- 
physical " ghost contribution" . More recently, Gould and 
Dobson [35] have made similar observations of " ghost in- 
teractions" in the context of the exact-exchange (EXX) 
functional with fractional spin densities, and used ensem- 
ble definitions to propose an improved, linearized EXX 
functional. 

Here, we offer an ensemble generalization of all en- 
ergy terms of an arbitrary density functional, to systems 
with a fractional N. Using the simplest functional of all - 
the local density approximation (LDA) - on example sys- 
tems, we find that this generalization greatly reduces the 
problem of the energy curve convexity, significantly re- 
stores the IP theorem, and concomitantly introduces an 
appropriate derivative discontinuity into the exchange- 
correlation potential in a natural manner. All this is 
achieved while neither introducing empiricism nor chang- 
ing the underlying functional form. 

Our considerations start with the ground state of a 
zero-temperature interacting-electron system with a frac- 
tional N, described by an ensemble (mixed) state of 
the form A = (1 - a)\^ No )(^ Na \ + a\^ No+1 )(^ No+1 \, 
where I'Jjvo+p) i s a (pure) many-electron ground state 
with No +p electrons and p is or 1 [44] . The electron 
density is then obtained using the many-electron density 
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operator, h(r) = Ei<K^ — as 

n(r) = Tr{An} = (1 - a)n {r) +an 1 (r). (2) 

no(f) and ni(f) are the densities of the interacting sys- 
tems with N and Nq + 1 electrons, respectively. As a 
result, the total energy E is obtained as in Eq. (1). 

In the Kohn-Sham (KS) formulation of DFT, the 
interacting-electron system is mapped into one KS sys- 
tem of non-interacting electrons with a fractional number 
of particles, N. Therefore, its ground state must also be 
an ensemble state, given by Aks = (1 — a )\^N^)^N^\ + 
a|$^ o ) +1 )($^ ) ) +1 |, where l^+p) are pure KS ground 
states, with No + p electrons, respectively [2] [44] [45] . 
Each of the pure ground states is described as a Slater de- 
terminant of single-electron orbitals {f\ a ^}- In contrast 
to the quantities \^n + p ) and n p of the many-electron 
interacting system, all quantities of the KS ensemble are 
a-dependent, a fact we emphasize via the superscript ( a >. 
Hence, in addition to the explicit dependence of Aks 
on a, there also exists an implicit dependence through 

Similarly to Eq. (2) for the interacting-electron sys- 
tem, the density in the KS system is obtained as 
ng(r) = Tr{k KS n} = (1 - a)p Q a \r) + ap ( °\r) = 
ESirfVT, where p p (f) := <#$ +p |n|$g +p ) = 
Et°i +P k W (0| 2 ,and 
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i^N 
i = N + l 
i > N Q + 1 



are the occupation numbers of the KS levels. Thus, the 
original ensemble is mapped into a KS system with a 
fractionally occupied frontier orbital. While n^j(r) is 
required to equal n(r) by construction, we stress that 
pp (r) need not equal n p (r). In other words, the fic- 
titious densities p p a \r), obtained by considering the 
frontier orbital as completely empty or completely full, 
need not correspond to the true electron densities of the 
integer-electron interacting systems, n p (r). Moreover, 
because no(r), n\{r) and n(f) can all be obtained inde- 
pendently from each other by considering systems with 
different (integer or fractional) electron numbers, Eq. (2) 
can be viewed as a point-wise linearity criterion for the 
density, complementing Eq. (1). 

We now examine the ensemble properties of the kinetic, 
external, and Coulomb energies of the KS system, associ- 
ated (in atomic units) with the operators T = — \ ^ i Vf , 

V n = v n (r)h(f), and W = \ £\ YLfri Vi - ^l" 1 , respec- 
tively. 

By inspection, the kinetic energy functional, 



TksH = (1 - a)<*$|f + a(^> +1 \f |*$ +1 > = 
Y^iLi 9i{f^\ ~ \^ 2 \ L P^)^ an( i the external energy 
functional, V n [n] = J v n (r)n^g( , r)d 3 r, retain their usual 
KS form and are explicitly linear in a. By definition [3], 



the Coulomb functional W = Ti{A KS W} = W H + W X is 
comprised of a Hartree (H) and an exchange (x) term. 
Performing the Tr operation and sorting out the Hartree 
and exchange contributions, we can express the ensemble 
terms Wh and W x by means of the standard, pure-state 
definitions of the Hartree and exact exchange (EXX) 
functional (see Supplementary Material for details). 
We obtain: 



W H = (l- a )E H [p ( Q a) }+aE H [p'; 
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where as usual 



{l-a)E x [p Q a) ]+aE x [p\ 
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and 



E x [n] 
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// 



d 3 rd 3 r 



,<P l i(r')<p* j (i f )<pi(T f )<p j (r') 



r — r 



(7) 

Because En[n] is not linear in n, it immediately follows 
that the required Wh of Eq. (4) is not obtained by in- 
serting the fractional-electron density n^j into Eq. (6). 
A similar statement is true for E x [n] and W x [35] . There- 
fore, the Hartree and EXX junctionals do not retain their 
usual form for ensemble states. Instead, Wh = Ejj[n] + 



(3) AE eH [^ +1 ;a] and W x = E x [n] - AE eH W 



where 



JVo+l' 



AE, 



eH 



-a(l 



j 3 ri 3 ; . >&i( f )| 2 |A ) + l^)| a 



r — r 



(8) 

is the ensemble (e) correction. 

Note that for a = or 1, Wh and W x reduce to their 
usual forms (6) and (7). Thus, introduction of the term 
AE e H should not affect the total energies of systems with 
an integer number of electrons. In addition, the total en- 
ergy obtained for EXX calculations with no correlation 
should not be affected cither, as AE e n appears with op- 
posite signs in Wh and W x [46]. However, the Hartree 
expression is usually complemented by an approximate 
xc-functional, E xc [n], that is not the EXX. Error can- 
celation is then not expected and, as shown below, not 
obtained. Trivially, an arbitrary £7 xc [n] is not linear in 
n, but it can still be made explicitly linear in a, in the 
same spirit as Eqs. (4), (5) above, yielding: 



E exc [n] = (1 - a)E xc [p a) ] + aE xc [p^']. 
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(9) 



(see additional argumentation in the Supplementary Ma- 
terial). For the special case of the LSDA, we refer to its 
ensemble generalized form, using Eq. (9), as eLSDA. 

Importantly, the ensemble expressions Wh (Eq. (4)) 
and E exc (Eq. (9)) no longer depend explicitly on the 



3 



density n, even for underlying functionals that are explic- 
itly density-dependent for pure states, such as the LSDA. 
Ultimately, they depend on the KS orbitals (themselves 
a functional of n) via (r) , as well as on a itself. This 
affects the KS potential, Vks- To remain within the 
KS framework, it must now be evaluated using the opti- 
mized effective potential (OEP) procedure, appropriate 
for implicitly density-dependent functionals [19, 36-38]. 
A complete derivation of vks is provided in the Sup- 
plementary Material. One unusual aspect of it, which 
we stress here, is that the explicit dependence of Wh 
and E exc on a contributes a spatially-uniform but in- 
dependent term to vks, given by 
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E xc y?>\-E a 



where v xc = 8E XC / 5n is the usual xc-potential. This term 
involves the highest (possibly partially) occupied orbital, 

Viv Vi' anc ^ does not vanish even when N is an inte- 
ger, despite the fact that for integer values the conven- 
tional and ensemble-generalized energy expressions are 
identical. Such a constant term, although allowed by the 
Hohenberg-Kohn theorem [39], is usually deemed unim- 
portant because it does not affect the density or the to- 
tal energy. However, it does shift the KS eigenvalues, a 
fact we show below to be crucial. Thus, all calculations 
now conceptually involve orbital-dependent functionals, 
although for integer N the term ti(°) can be easily evalu- 
ated without performing the computationally demanding 
OEP calculation. 

To illustrate the proposed generalization and its im- 
plications, we apply the eLSDA functional to the H2 
molecule and the C atom - the latter being frequently 
used to demonstrate the problem of convexity in the en- 
ergy versus fractional charge curve [9, 14, 16, 18]. All 
calculations were performed within the Krieger-Li-Iafrate 
(KLI) approximation to the OEP solution [38] , using the 
all-electron, real-space DARSEC code [40]. Numerical 
details are given in the Supplementary Material. 

The total energy for the above two systems, as a func- 
tion of the net charge, q, are given in Fig. 1, with q 
ranging from -2 (doubly-ionized system) to (neutral 
system) . The LSDA energy curves are, as expected, con- 
vex [9, 14, 16, 18]. The curve for the eLSDA is, however, 
almost piecewise linear, being slightly concave. The su- 
periority of eLSDA over LSDA in piecewise-linearity is 
quantified by measuring the curvature of E(q), given in 
Table I [4 7]. Clearly, the strong reduction in the devia- 
tion from piecewise-linearity is a significant advantage of 
the ensemble approach. This deviation is not fully elim- 
inated because, while the eLSDA functional is explicitly 
linear in a by construction, it may still be implicitly non- 
linear as the KS orbitals themselves are also a-dependent. 
Comparison of the eLSDA results to the EXX ones, also 
given in the Figures and Table, show that the piecewise- 
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FIG. 1: Energy of the H2 molecule (top) and of the C atom 
(bottom) as a function of fractional charge q, for various func- 
tionals. EXX results for H2 have been shifted upwards by 0.4 
Ry, for clarity 
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TABLE I: Curvature of E(q), in Ry, obtained for H2 and C 
with various functionals 



linearity of eLSDA is comparable to that of EXX. An 
obvious advantage of eLSDA, however, is the treatment 
of correlation. 

An additional perspective on the improvement afforded 
by eLSDA is provided by considering the extent to which 
the density linearity criterion, Eq. (2), is obeyed. Wc 
consider D(f) := n(f) — (1 — a)no(r) — an\{r), i.e., the 
difference between the left hand and right hand sides of 
Eq. (2), which should equal at all f for the exact func- 
tional. A plot of D(r) at q = —0.5 for H 2 , as obtained 
with LSDA and eLSDA, is presented in Fig. 2. Clearly, 
the spatial profile of D{r) is smoother with eLSDA and 
its average numerical value much smaller. Specifically, 
Q(q) := / D 2 (r)d 3 r, which is the variance of D{r) per a 
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FIG. 2: Deviation from piecewise-linearity in the density, 
D(r), obtained for the H2 molecule for q = —0.5 using the 
(a) LSDA, (b) eLSDA 



given q, is Q ~ 1CT 4 Bohr- 3 with LSDA. With eLSDA, 
however, it is lower by two orders of magnitude for 
q = — 1...0 and essentially zero for q = —2... — 1. 

The great improvement in the piecewise-linearity of the 
energy curve (Fig. (1)) is directly manifested in the de- 
gree to which the IP theorem is satisfied. This is illus- 
trated in Fig. 3. The figure shows the highest (possi- 
bly partially) occupied orbital, £^ D , the energy deriva- 
tive dE/dq, and the negative of the ionization energy, 
—I (computed from total energy differences obtained at 
integer q values), as calculated for H2 as a function of 
q with both LSDA and eLSDA. Janak's theorem [41], 
which equates between and dE/dq for any approxi- 
mate functional, is indeed closely obeyed by both approx- 
imations, confirming the numerical accuracy of our cal- 
culations. But because eLSDA is much more piecewise- 
linear, EhoW) calculated with it is much more piecewise- 
constant as a function of q (as it should be for the ex- 
act functional). Furthermore, Eho coincides much more 
closely with —I when approaching an integer q from be- 
low, in agreement with the IP theorem [8, 13]. Thus, at 
integer occupations the IP theorem is much more closely 
obeyed - a direct consequence of the above-mentioned 
constant potential v^°\ even though the the total energy 
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FIG. 3: Frontier orbital energy, Eho, energy derivative dE/dq, 
as a function of q, as well as the negative of the ionization en- 
ergy —J, calculated for H2 with the LSDA and eLSDA func- 
tional 
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'For H^~, no experimental value for I2 exists. Instead, it was 
obtained from EXX calculations, which yield an exact result for 
this system. 

TABLE II: The highest occupied orbital energy —Eho, com- 
pared to the ionization energy I, for the neutral H2 and C, as 
well as the cation fundamental gap, deduced from the discon- 
tinuity of Eho at q = — 1, compared to the difference between 
the second and first ionization energies of the neutral system, 
obtained with various functionals. All energies are given in 
Ry. A's correspond to the relative error between the two 
values positioned immediately above them. 



has not changed. 

The satisfaction of the IP theorem is closely related to 
another fundamental property of the exact xc-functional: 
as the number of electrons crosses an integer, the xc- 
potential may "jump" by a constant, usually known 
as the derivative discontinuity (DD) [8]. The conven- 
tional wisdom on explicit density functionals (includ- 
ing LSDA) is that they do not possess this discontinu- 
ity. Instead they "average over" it and consequently 
Eho is found to deviate from — / by about half of the 
missing DD [43]. Recently, Stein et al. [12] have shown 
that the missing DD and the deviation from piecewise- 
linearity are doppelgangers: a significant increase in the 
degree of piecewise-linearity must be accompanied by 
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the appearance of a discontinuity in the xc-potential. It 
emerges from the above-discussed constant term, v(°\ 
which arises from the orbital-dependent treatment of the 
KS potential, depends on the highest (possibly par- 
tially) occupied orbital and is therefore different if one 
approaches an integer N from the left or from the right. 
Therefore, the DD of explicit density functionals arises 
naturally, without invoking any empiricism, and is calcu- 
lated at a negligible computational cost. This is readily 
observed in Fig. 3 and Table II: the fundamental gap of 
the ion Hj, as deduced from the difference of e^o to the 
left and right of q = —1, is much larger with eLSDA 
(1.320 Ry) than with LSDA (0.426 Ry), and corresponds 
much more closely to the result obtained from total en- 
ergy differences (1.489 Ry, solid black line in the figure). 
Similar observations apply to the C atom (see Table). 
Thus, our ensemble-based approach automatically iden- 
tifies and restores the missing derivative discontinuity, 
appropriate for the underlying functional. 

In conclusion, we presented a generalization of the 
Hartree, exchange and correlation terms of an arbitrary 
density functional to systems with a fractional electron 
number, based on the ensemble form of DFT. Using the 
local density approximation on H2 and C as illustra- 



tive examples, we showed that this generalization signifi- 
cantly reduces the deviation from piecewise linearity and 
generates the appropriate derivative discontinuity, with- 
out introducing empiricism and with no changes to the 
underlying functional form. With this generalization, the 
total energy at integer electron numbers remains intact, 
but the eigen-energies change and the ionization potential 
theorem is much more closely obeyed. This shows that 
problems that have plagued simple approximate density 
functionals for many years can be very strongly miti- 
gated by rigorous employment of ensemble DFT within 
the optimized effective potential approach, without any 
further functional development. We expect this proposed 
generalization to be equally useful for more advanced ap- 
proximate functionals, as well as for more complex sys- 
tems, allowing for improvement in spectroscopic proper- 
ties without any compromise on energetics. 
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In the following, reference to equations from the main 
text are made by a single number (e.g. Eq. (1)), while 
equations from this document include the section num- 
ber, as well (e.g. Eq. (1.1)). 

I. DERIVATION OF EQS. (4), (5), (8) 

To obtain Eqs. (4) and (5), consider the application of 
the Coulomb operator 



N N 



i=l j = l 11 -? 1 



(1.1) 



to a pure state |$) of the KS system, where 



¥>i(rl) tpi(f2) •■• ipiirw) 
^jv(rl) (pN(ri) ■■■ ¥>iv(riv) 



(1.2) 



is a Slater determinant of the one-electron KS orbitals 
{tpi}. The symbol ^ does not accompany here the KS 
orbitals and other derived quantities, as opposed to the 
main text, for clarity of presentation. 

The quantity W := ($\W\$) is given by the well- 
known expression 



N N 



w 
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N N 

»=1 j=l J J 



,3 .3 ,<Pi(f)<Pi(f)<Pj( r ')<Pj(. r ') 

a rd r — 



-3 .3 , ( Pi( r ') ( Pj(r)¥>i(r) l Pj( r ') 
d rd r ^ 



(1.3) 



In the first term above we recognize the pure-state 
density, defined as p — J2iLi I | 2 3 an d partition the 
Coulomb energy as 



W = E H [p\+E x [p], 



where 



Eh[p] = \ If d 



'3 rd 3 r/ P(^P( r ') 



\r — r 



(1.4) 



(1.5) 



is the Hartree energy, and 



N N 



(1.6) 

is the exchange energy (where ipi themselves are func- 
tional of the density). 

For an ensemble state, W is obtained using the Tr 
procedure: W = Tt{AksW}, where 

A^s = (1- a)\$ No )(<f> No \ +a\®N +i)(<S>N +i\, (1.7) 
and 

W= (I - a)($ Na \W\$ No ) + a(<S> Na+1 \W\$ No+1 ). (1.8) 

The densities corresponding to the states |$jv ) and 
|$at +i) are po and pi, respectively (see definitions in 
the main text), and therefore 

W = (l-a)E H \M+^ E H\pi] + 0--oi)E x \po]+aE x \p 1 ]. 

(1.9) 

The first two terms are the Hartree energy of the ensem- 
ble state, denoted by Wr- The last two terms are the 
exchange energy of the ensemble state, denoted by W x . 
These four terms appear in Eqs. (4) and (5). 

To express Wh in terms of E H [n] and AE eH , we insert 
the definition n = (1 — a)po + api (Eq. (2)) into Eq. (1.5) 
to obtain 



E H [n] = (1 - a) 2 E H [p a ] + a 2 E H [ Pl ] 

+ a{\ - a) II d s rd'r 
= (1 - a)E H [p ] + aE H [pi] 



., l(1 _ fl)| jjd W P -f^l 



E H [po] - E H [pi 



(1.10) 
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Combining Eqs. (1.5) and (1. 10) while using the fact 
that pi - po = \<Pn +i\ 2 , yields 



E H [n] = (1 - a)E H [p ] + aE H [ Pl ] 

i r r . . 

- -a(l -a) d 3 rd 



/\|2 
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Thus, we arrive at the relation Wh = E H [n] + 
AE eH [ip No+1 ;a\. 

To express W x in terms of E x [n] and AE e n, we rear- 
range Eq. (5) as 

W x = E x [p ] + a(E x [ Pl ] - E x [p }). (1.12) 

Using Eq. (1.6) one then obtains 

E x [ P i] - E x [p Q ] = 



_ 1 ^ || ^^a^ ^o+i^O^^^o+i^^^) 



r — r 



(1.13) 



From Eq. (3) 9n +i = a; considering Eq. (1.12), while 
substituting Eq. (1.6) for its first term and Eq. (1.13) for 
the second term, we obtain 
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Finally, to achieve the form of Eq. (7), we replace the 

sum by Y^jLi 9i- This change introduces an addi- 

tional term to the double summation, which should be 
subtracted for maintaining the equality. This manipula- 
tion leads to the expression 



= -*£!>&• ii d " rd 

%=\ j=i j j 



\r — r 

d 3 rd 3 r , ^No+lir'^^No+lir)? 

\r — r'\ 

(1.15) 



which is equivalent to the form W x — E x [n] 
AE eH [(p No+1 ;a\. 



II. RATIONALE OF EQ. (9) 

In the following we explain in more detail the rationale 
behind the ensemble generalization of an approximate xc- 
functional to a form which is explicitly linear in a. 

The approximate xc-functional can be (and usually is) 
presented as a sum of an approximate exchange func- 
tional and an approximate correlation functional. Be- 
cause the exact exchange functional is explicitly linear in 



a, as shown in Eq. (5), it is reasonable to require that 
the approximate exchange functional exhibit the same 
property. 

The correlation functional for a pure state can be for- 
mally expressed, without any approximation, as (see e.g. 
Eq.(1.68) in Ref. [1]) 



E c = (*|T + W\V) - ($|T + W\§), 



(11.1) 



where T and W are the kinetic and the Coulomb opera- 
tors, respectively, \& is the wavefunction of the interacting 
system, and $ is the wavefunction of the KS system. 

For ensembles, the states |<J/) and |<&) are represented 
by the operators A and Aks, respectively (see definitions 
in the main text). The correlation energy then becomes 

E c = (l- a)(* No \f + W\V No ) + a(V Na+1 \f + W\* No+1 ) 

-(1 - a)<*g|T + W\*%) + a<*$ +1 |T + W\*% +1 ), 

(11.2) 

being explicitly linear in a. Thus, it is reasonable to 
require that the approximate correlation functional will 
also be explicitly linear in a. 



III. DERIVATION OF THE KS POTENTIAL 

In our ensemble approach, the KS potential, vks, can- 
not be obtained by a straightforward application of a 
functional derivative, because parts of the energy func- 
tional E[n] - T KS + V n [n] + E H [n] + &E eH [(p$ +1 ; a] + 

Eexc [Po^ , Pi"' 1 ; ol] depend explicitly on the orbitals and 
a, rather than on n. 

Let us denote vks = V n+ V H + v eH + v exc , where v H = 
8E H /8n, v eH = 5(AE eH )/5n and v exc = SE exc /Sn. The 
two last terms of vks can be presented as 



v T (r) 



SE T 
5n 



dE- 



t \ 5a ( 5E T \ / ttt -, \ 

srj.srnsw).' (IIL1) 



where T stands for either the eH term or exc term. 
For the ensemble of two pure, non-degenerate states dis- 
cussed throughout, a[n] — N-Root(N) and N = J nd 3 r. 
We then find 5a/ Sn = 1. The second term of the RHS 
of Eq. (III.l) is denoted by and can be expressed as 



oo 

i=i 



SEh 



5(fi(r') 



+ c.c, 



(III.2) 

where a is constant and therefore treated as a parameter. 
Because dEx/Sipi can be obtained from Eqs. (8) and (9), 
can be found using the OEP procedure [2], as is 
suitable for an orbital-dependent functional. 

The first term of the RHS of Eq. (III.l) is denoted 

by . This term is somewhat unusual, as most func- 
tionals do not contain the quantity a explicitly, is 
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a-dependent, but constant in space. As the expressions 
for AE e H and E exc do not depend explicitly on the den- 
sity n, it is not possible to take the derivative (dEx / da) n 
directly. However, the quantity (8Et / da){ Vi } is acces- 
sible. To relate between the two quantities, we express 
Et as a functional of a and n, where the latter is itself a 
functional of a and {ifi}: Et — Et [a, n[{tpi}, a]]. Then, 
we obtain 



dE T 
da 



0E T 
da 



dn(r) 
da 



(III.3) 

On the right-hand side of this expression we recognize the 



first term to be and v T L ' to be the first multiplicand of 

the second term. From Eq. (2), (dn/da){ Vi y = l^jvo+il 2 - 
As a result, 



,(!) 



dE 7 
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- / i^ViWiH 15 ^^. (IH.4) 

{•Pi} J 

We remind that in the OEP formalism for the highest 
occupied orbital 

:=/kg +1 (r)| 2 4 1) (^ 3 r = 
= / \^l +i m 2 UT{^r=-.u T , Na+l , (III.5) 



where 
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(see [2], Sec.ll and [3]). An analytical form for VjP can 
then be obtained from Eqs. (8) and (9). Using this rela- 
tion, for the eH functional we find: 
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and for the errc-functional: 
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-d 3 r, (III.7) 
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(III. 



During the self-consistent numerical solution of the KS 
equations with the proposed functional, the spatial con- 
stant ' can be omitted, as the addition of a constant 
to the potential does not affect the eigenfunctions, the 
density, or the total energy E. However, has to be 
taken into account when addressing KS eigenenergies, in 
particular when comparing the energy of the frontier or- 
bital to the derivative dE/dq (see definition in the main 
text), following the Janak and the IP theorems [4, 5], and 



when calculating the energy gap as varies with No 
and thus does not cancel out. 

The equations above can be generalized to the spin- 
polarized well. Because in the spin-polarized 
version of DFT there exist two potentials, Vks,<t, where 
(j =t or |, there also exist two sets of orbitals, {f[ a J}, 
two densities, n CT , and also two statistical weights a a . As 
a result, Eq. (III. I) is generalized to be 



v T,a(r) 
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5n a 



'dE T \ f SE T 
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(III.9) 

where r refers to the other spin channel than a. Other 
equations of this section can be generalized accordingly; 
In particular, Eq.(III.4) reads: 



(0) _ / dEj 



VT >°-\da a , 

{ Vi , T },a T 

(111.10) 

We stress, however, that in all formalism presented in 
the main text and here only the a of one spin channel is 
allowed to be fractional. 



IV. NUMERICAL DETAILS ON CALCULATION 
OF THE SYSTEMS H 2 AND C 

The bond length of H 2 was found by relaxation to be 
1.45 Bohr for LSDA (and therefore, by definition, also for 
eLSDA), and 1.39 Bohr for the EXX functional. This re- 
sult is in close correspondence with previous calculations 
(see Ref. [6] and references therein) . The bond length was 
kept unchanged when varying the number of electrons in 
the system. 

Calculations with the orbital-dependent functionals, 
cLSDA and EXX, were performed using the Krieger-Li- 
Iafratc (KLI) approximation [7] , which is less demanding 
computationally. 

In the calculations performed for the C atom, it was 
assumed that the neutral C has a spin S z = 1, the ion 
C + has the spin of S z = |, and for the ion C ++ , S z = 0. 
Therefore, varying the number of electrons in the system 
was performed solely in the spin-up channel. In addition, 
the axial quantum number L z was restricted to be for 
q = —2... — I, and increased linearly with q for q = — 1 ...0, 
to obtain L z = 1 for the neutral C atom. Calculations 
with other values of L z were checked as well: the total 
energy they produced differed from the reported values 
by less than 4 mRy. These restrictions for the C system 
assured that the calculation is performed with an ensem- 
ble of two states, which was the one considered in the 
main text. 
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